Variational ground states of 2D antiferromagnets in the valence bond basis 
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We study a variational wave function for the ground state of the two-dimensional S = 1/2 Heisen- 
berg antiferromagnet in the valence bond basis. The expansion coefficients are products of ampli- 
tudes h(x, y) for valence bonds connecting spins separated by (x, y) lattice spacings. In contrast 
to previous studies, in which a functional form for h(x, y) was assumed, we here optimize all the 
amplitudes for lattices with up to 32 x 32 spins. We use two different schemes for optimizing the 
amplitudes; a Newton/conjugate-gradient method and a stochastic method which requires only the 
signs of the first derivatives of the energy. The latter method performs significantly better. The 
energy for large systems deviates by only w 0.06% from its exact value (calculated using unbiased 
quantum Monte Carlo simulations). The spin correlations are also well reproduced, falling w 2% 
below the exact ones at long distances. The amplitudes h(r) for valence bonds of long length r 
decay as r -3 . We also discuss some results for small frustrated lattices. 

PACS numbers: 75.10.Jm, 75.10.Nr, 75.40.Mg, 75.40.Cx 



I. INTRODUCTION 

The valence bond (VB) basis for singlet states of quan- 
tum spin systems was first discussed by Rumer— and 
Pauling 2 - in the 1930s and was shortly thereafter applied 
in Heisenberg spin chain calculations by Hulthen. 3 For a 
systems of N S = 1/2 spins the overcomplete and non- 
orthogonal VB basis consists of states that are products 
of N/2 spin pairs forming singlets. In the most general 
case the members of a singlet can be separated by an ar- 
bitrary distance, but it is often convenient to consider a 
restricted set of configurations with only bonds connect- 
ing two different groups of sites, e.g., the two sublattices 
of a bipartite lattice. Such a restricted VB basis is still 
overcomplete and any singlet state can be expanded in 
it. Any restriction on the maximum length of the bonds 
will render the basis incomplete, however. 

After Anderson's proposal in 1987 of a resonating- 
valence-bond (RVB) state 4 as a natural starting 
point for understanding high-temperature supercon- 
ductivity in the cuprates^ variational VB states 
were investigated for both doped and undoped 
antiferromagnets^i^'iSiii^ii 3 . The RVB spin liquid 
mechanism is based on states dominated by short va- 
lence bonds, which in the extreme case have been argued 
to correspond closely to the quantum dimcr model. 14 
It was early on established, however, that the ground 
state of the two-dimensional (2D) Heisenberg model 
with nearest-neighbor interactions, which has a Nccl or- 
dered ground state and describes very well the undoped 



cuprates 



15,16 



actually requires an algebraic, not expo- 



nential, decay of the bond-length probability. 6 

Some attempts were made to use the VB basis as a 
framework for numerical calculations ) 10 ' 17 : 18 ' 19 : 20 ' 21 but, 
with very few exceptions^ these efforts were not pur- 
sued further in large scale quantum Monte Carlo cal- 
culations (QMC). However, it was recently pointed out 
that there are previously unnoticed advantages of carry- 
ing out ground state projector QMC calculations in the 



VB basis, including a natural way to access excitations in 
the triplet secto n 23 ' 24 ' 25 Such a scheme has already been 
applied to 2D and 3D models with valence bond solid 
ground states i 26 ' 27 

It may also be worthwhile to pursue further variational 
schemes in the VB basis, especially exploring possibilities 
to study frustrated systems this way. In this paper we re- 
port a bench-mark variational calculation going beyond 
previous VB variational studies 6 - of the 2D Heisenberg 
model. Instead of assuming a functional form for the 
bond-length amplitudes and optimizing a few parame- 
ters, we optimize all individual amplitudes in order to 
definitely establish the properties of this kind of wave 
function and its ability to reproduce the ground state of 
the 2D Heisenberg model. We also report some prelimi- 
nary studies of a frustrated system. 

For the standard 2D Heisenberg model with nearest- 
neighbor coupling J, the energy of our best optimized 
wave function deviates by only AE/J w 0.06% from the 
exact ground state energy for system with up to 32 x 32 
spins. The size dependence of AE shows that this ac- 
curacy should persist in the thermodynamic limit. The 
error is only half that of the best wave functions found 
in the previous variational QMC study— The spin-spin 
correlations are also remarkably well reproduced; they 
are approximately 2% smaller than the exact values at 
long distances, corresponding to an w 1% underestima- 
tion of the sublattice magnetization. We find that the 
asymptotic form of the amplitudes for bonds of length r 



is h(r) 



which has also been found recently in a 



mean- field calculation. 28 We also compare the variational 
wave function with the exact ground state in the case of 
a 4 x 4 lattice, and find that the overlap is w 0.9998. Ex- 
tending the 4x4 calculation to a frustrated system, the 
Ji— J2 model, we find that the quality of the amplitude- 
product wave function deteriorates as the frustration is 
increased but the overlap remains above 0.996 even for 
J2/J1 as high as 0.4. 

In Sec. II we define the variational VB wave function. 
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Technical details of the QMC based optimization meth- 
ods that we have used to minimize the energy are pre- 
sented in Sec. III. We discuss both a standard Newton 
method and a stochastic scheme, which requires only the 
first derivatives of the energy. The latter method per- 
forms significantly better for large lattices. Results for 
the energy and spin correlations in the standard non- 
frustrated Heisenberg model are discussed in Sec. IV. 
Results for energies and overlaps for the 4x4 frustrated 
lattice are discussed in Sec. V. In Sec. VI we conclude 
with a brief summary and discussion of the methods and 
results. 



II. MODEL AND WAVE FUNCTION 



We study the standard 5=1/2 Heisenberg model; 



H — J Si ■ Sj, 

(id) 



(1) 



where denotes nearest-neighbor sites on a 2D square 
lattice and J > 0. The basic properties of this model have 
been known for a long time^ and ground state param- 
eters such as the sublattice magnetization, the energy, 
and the spin stiffness have been extracted to high pre- 
cision in many QMC studies! 21 ' 29 ! 30 Here our aim is to 
investigate how well a simple variational wave function 
can reproduce the true ground state. 

The general form of a VB wave function for N S = 1/2 
spins is 



l*> - EM *' 6 ?)' 

k 



( a AT/2i b%/ 2 )) 



(2) 



where (a*, b^) represents a singlet formed by the spins at 
sites a and b in VB configuration k; 



(a,b) 



V2 



(lalb — lalb)- 



(3) 



The notation |V&) has been introduced in @ for conve- 
nience. In the most general case, the VB configurations 
Vfc include all the possible pairings of the N spins into 
N/2 valence bonds. A more restricted, but still massively 
overcomplete basis is obtained by first dividing the sites 
into two groups, A and B, which in the case of a bipartite 
lattice naturally correspond to the two sublattices. Here 
we will use such a restricted VB basis and always (also 
when considering the frustrated case later on) take A and 
B to refer to the sublattices of the square lattice. We fix 
the "direction" of the singlet in © by always taking the 
first index in (a, b) from A and the second one from B. 
With this convention one can show that all the expansion 
coefficients [where k = 1, . . . , (JV/2)!] in Eq. |2]) can 
be taken positive. This corresponds to the Marshall sign 



rule for a non-frustrated system in the basis of eigenstates 
of the operators.^ - 

We consider expansion coefficients of the amplitude- 
product form previously introduced and studied by 
Liang, Doucot and Anderson^ 



N/2 N/2 

f k = l[h(alb^ = l[h(xly^, 



(4) 



(=1 



where X4 and yi are the x and y separations (number 
of lattice constants) between sites at,bi, which are con- 
nected by valence bond i. Considering the lattice sym- 
metries (we use periodic boundary conditions), there are 
hence w AT/16 independent amplitudes h(x,y) to opti- 
mize. In the previous study^ only a few short-length 
amplitudes were optimized and beyond these an asymp- 
totic form depending only on the length r of the bond 
was assumed^ With a power-law form, h(x, y) ~ r~ p , it 
was found that long-range Neel order requires p < 5. The 
best variational energy was obtained with p = 4, giving 
a deviation AE/J rj 0.0008 (or w 0.12%) from the exact 
ground state energy (the thermodynamic-limit value of 
which is^ 9 - w 0.69944 per site). 

We here optimize all h(x, y) using two different meth- 
ods: A standard Newton method (combined with a con- 
jugate gradient method — the Fletcher Reeves method 32 ), 
and a stochastic method that we have developed which 
requires only the signs of the first derivatives. 



III. OPTIMIZATION METHODS 

We now discuss the technical details of these calcula- 
tions. To optimize the energy using a Monte Carlo based 
scheme, we write its expectation value as 



E = (y\Hm = 



(Vk\H\Vi) 
Wk\W) 



EufkMVklVt) 



(5) 



where on the right-hand side we have taken into account 
the fact that we do not normalize the wave-function co- 
efficients, i.e., the amplitudes in Eq. (|3|) are only de- 
termined up to an over-all factor. The overlap (Vfe | VJ) 
between two VB states is related to the loops form- 
ing when the two bond configurations are superimposed; 
(V k \Vi) = 2 N '- N / 2 , where N t is the number of loops^SS 
Matrix elements of Si • Sj are also easy to evaluate in 
terms of the loop structure. If i and j belong to the 
same loop, then (V4|S, • S rf | V*>/<V* | V?> = ±3/4 (+ for i,j 
on the same sublattice and — else), and else the matrix 
element vanishes Recently more complicated matrix 
elements have also been related to the loop structure! 2 ^ 
For a given set of amplitudes h(x,y), we evaluate the 
energy using the Metropolis Monte Carlo algorithm de- 
scribed in Ref. @. An elementary update of the bond 
configuration amounts to choosing two next-nearest- 
neighbor sites a, c (or in principle any a, c in the same 
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sublattice, but the acceptance rate decreases with in- 
creasing distance between the sites), and reconfiguring 
the bonds (a,b),(c,d) to which they are connected, ac- 
cording to (a,b)(c,d) — > (a, d)(c, b) (where the order of 
the labels here correspond to both sites a and c being in 
sublattice A). The Metropolis acceptance probability P 
for such an update is very easy to calculate in terms of 
amplitude ratios and the change in the number of loops, 
ANi, in the overlap graph; 



P 



HXgd, yad)h(Xcb, Vcb) 2 AjV ; l 

h(x ab ,y ab )h(x cd ,y cd ) 



(6) 



A. Newton Conjugate Gradient Method 

For the optimization we also need derivatives of the en- 
ergy with respect to the amplitudes. The Newton method 
requires first and second derivatives. Moving in a certain 
direction g in amplitude space, the amplitude vector h n 
is updated from iteration n to n + 1 according to 



hn+i — h„ — 



^(hn) 

E"(h n ) 



g. 



(7) 



where E' g and E 1 ' are the first and second derivatives 
of the energy along the g direction. They are calcu- 
lated from the derivatives with respect to the amplitudes 
h(x,y), which are evaluated during the sampling of VB 
configurations. Writing the energy expectation value as 



(E) 



J2 p w p e p 



W p = ~[[h(x,y) n ° 



(8) 



where A is the Hessian matrix. In practice, we have 
found that the number of line optimizations required for 
energy convergence is of the same order as total number 
of different amplitudes. Since the optimization is based 
on quantities obtained using a stochastic scheme, the fi- 
nal results of course have statistical errors. 

We have oc A^ 2 different second derivatives and hence 
computing all of them requires a significant computa- 
tional effort. Their statistical fluctuations are also rel- 
atively large for large lattices. We have also used a 
minimization method requiring only the first derivatives, 
where the amplitudes are iterated according to 



E' g (h n )h n ^ ~ E' g (h n ^)h n 



(13) 



This amplitude update, which gives the exact minimum 
in the quadratic regime, works well when h(x, y) are get- 
ting close to their optimum values. However, we still need 
the second derivatives until we get very close to the opti- 
mum and in practice the advantage of using (|13p instead 
of ([7]) at the final stages does not appear to be significant. 

We here note that a method recently proposed to re- 
duce the statistical fluctuations of the Hessian in varia- 
tional QMC optimizations of electronic wave functions 33 
is not applicable here. The proposal was to symmetrize 
the Hessian and write it solely in terms of covariances, 
by adding terms which are zero on average but reduce 
the statistical fluctuations of a finite sample. However, 
in our case the Hessian (jlOllip is already of this form 
and there is nothing more to do in this regard. 



where n xy is the total number of VBs of size (x, y) in the 
VB configurations Vk and Vi , the Monte Carlo estimator 
for the first derivative is 



d(E) _ / r hLF 



dh a \ h a 



(9) 



Here, to simplify the notation, we use a as a collective 
index for (x, y). The second derivatives — the elements of 
the Hessian matrix — are: 



1 



d 2 (E) _ 

dhl ~ hi 



\n\E) - (nl)(E) - (n a E) + (n a )(E) 



+2(n a ) 2 (E) -2{na)(naE)), 



(10) 



d 2 (E) 



1 



n a n b E) - (n a n b )(E) + 2(n a )(n b )(E) 

(11) 



dh a h b h a h b 

-(n a )(n b E) - (n b )(n a E) 



Since we have many amplitudes h(x,y) to optimize, we 
choose our optimization direction by the conjugate gradi- 
ent method. The first direction is the gradient direction. 
In subsequent steps we choose the direction conjugate to 
the former one, satisfying the relation 



ln+1 



A • h„ 



0. 



(12) 



B. Stochastic method 

We have developed a completely different optimization 
method which turns out to work much better than the 
Newton method. It is a stochastic scheme which only 
requires the signs of the first derivatives. We update 
each amplitude h a according to 



HK 



(14) 



where R is a random number in the range [0, 1) and 
sign(ir) = 1 if x > and —1 for x < 0. The parameter (3 
is gradually reduced, so that the amplitude changes be- 
come gradually smaller. If this "annealing" is performed 
slowly enough the scheme converges the system to the 
lowest energy. In fact, it turns out that the random fac- 
tor R is not even needed; the scheme converges even if 
R = 1. With R = 1, this kind of updating scheme has 
in fact been introduced previously in the context of neu- 
tral networks and is known as Manhattan learningi 34 ^ 5 
Note that even in this case there is still in principle some 
randomness in the method, because when the derivatives 
become small the signs of their QMC estimates can be 
wrong due to statistical fluctuations. Thus, there are 
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occasional adjustments of amplitudes in the wrong direc- 
tion. We have found that the random factor R speeds 
up the convergence, especially in the initial stages of op- 
timization, and so we normally include it. 

Our method is also related to what is known as 
stochastic optimization^-^ where the parameter vector 
is updated in a steepest-decent fashion according to the 
stochastically evaluated gradient; 

in(K +1 ) = HK)-0(^P} ■ (15) 

\ a / n 

This method has been used by Harju et al£& to opti- 
mize electrinic wave functions. We have tested it on the 
present problem but find that it performs significantly 
worse than our own variant of stochastic optimization. 
The reason appears to be that the fluctuations in the gra- 
dient can be very large, which can cause large detremen- 
tal jumps in the configuration space. In our scheme the 
step size is bounded and we avoid such problems. 

As in simulated annealing methods, a slow enough re- 
duction of (3 should give the optimum solution. In the 
present case, unlike in simulated annealing, the optimum 
reached should only be expected to be a local optimum, 
although the stochastic nature of the scheme does al- 
low for some more extensive exploration of the parame- 
ter space than with deterministic schemes. In stochastic 
optimization using the gradient, it has been argued that 
an annealing scheme of the form 



(16) 



should be used in order for the method to converge. We 
have found this scheme with a rj 3/4 to work well. 

It would be interesting to see whether this very sim- 
ple scheme could also be applied to optimize wave func- 
tions in electronic structure calculations — the time sav- 
ings from not having to calculate second derivatives are 
potentially very significant in problems with a large num- 
ber of parameters. 



IV. RESULTS 

We now discuss our results. In Fig. [1] we compare the 
ground state energy as calculated with the stochastic and 
Newton methods. We show the deviations from the cor- 
rect ground state energies (obtained using unbiased QMC 
calculations^) as a function of the lattice size L up to 
L = 32. Both optimization methods give very small en- 
ergy deviations for the 4x4 system — about 0.005% — but 
the error grows as the lattice size increases. We have cal- 
culated error bars by repeating the optimizations (from 
scratch) several times. It should be noted, however, that 
the fluctuations in this kind of nonlinear problem do not 
necessary have expectation value zero. Hence error bars 
calculated in the standard way in general only reflect 
partially the actual errors. 
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FIG. 1: (Color online) The deviation AE = E V3 , T {L) - 
-Eexact (L) of the variational ground state energy from the ex- 
act energy as a function of lattice size. Results obtained with 
both the Newton method and the stochastic opitimization 
scheme are shown. The error bars were obtained by carrying 
out several independent optimization runs for each L. 



For L > 10 the stochastic method delivers significantly 
lower energies, indicating that the Newton method has 
difficulties in locating the optimum exactly. We believe 
that this problem is due to the statistical errors of the 
second derivatives (which are much larger than those 
of the first derivatives). Convergence issues related to 
statistical fluctuations of the Hessian are well known in 
variational calculations for electronic systems^ Any sys- 
tematic shifts in the Newton method would of course be 
reduced by increasing the length of the simulation seg- 
ments used to calculate the energy and its derivatives in 
each iteration. Simulations sufficiently long to reach the 
same level of optimization as with the stochastic method 
do not appear to be practically feasible, however. In the 
stochastic scheme all statistical errors should decrease to 
zero as the cooling rate is reduced. We cannot, of course, 
guarantee that the results shown in Fig. [T] are completely 
optimal, but we have carried out the simulations at dif- 
ferent cooling rates in order to check the convergence. 
Based on these tests we do believe that the results are 
converged to their optimum values to within the error 
bars shown. 

For the larger lattices, the energy deviation in Fig. [1] 
is only AE/J w 0.0004, or w 0.06%, and is size in- 
dependent within statistical errors for L > 10. This 
should then be the accuracy in the thermodynamic limit. 
In Ref. only a few of the short-length amplitudes 
were optimized and a functional form — power-law or 
exponential — was used for the long-range behavior. The 
best power-law wave functions had energy deviations of 
AE/J w 0.0008; twice as large as we have obtained here 
with the fully optimized amplitudes. 

Having concluded that the stochastic method is the 
preferred optimization technique, we discuss only results 
for other quantities obtained this way. Fig. [2] shows the 
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L 

FIG. 2: (Color online) The staggered structure factor S(n, tt) 
versus lattice size, compared with unbiased QMC results. The 
inset shows the long-distance spin-spin correlation. Statistical 
errors are smaller than the symbols. 

size dependence of the staggered structure factor, 

S(-K,n) = Y l {-^ x+y C(x,y), (17) 

where C(x,y) is the correlation function, defined by 

C(xi - x 3 ,yi - y/) = (Si ■ Sj). (18) 

The inset shows the correlation function at the longest 
distance; (x,y) = (L/2, L/2),. We again compare with 
results from unbiased QMC calculations^ The structure 
factor of the variational ground state agrees very well 
with the exact result for these lattice sizes — the devi- 
ations are typically less than 0.5%. The long-distance 
correlations show deviations that increase slightly with 
L, going to « 2% below the true values for L > 10 
[which then should also be the asymptotic L — > oo error 
of S(it,it)]. The sublattice magnetization is the square- 
root of the long-distance correlation function; it is thus 
only 1% smaller than the exact value. 

Liang et al. did not conclusively settle the question 
of the asymptotic behavior of the amplitudes h(x, y) for 
bonds of long length r = (x 2 + y 2 ) 1 / 2 ^ The best vari- 
ational energy was obtained with an algebraic decay; 
h ~ r~ 4 . However, the energy is not very sensitive to 
the long-distance behavior and the values obtained for 
r -3 and r~ 2 were not substantially different at the level 
of statistical accuracy achieved. Even with an exponen- 
tial decay of the bond-length distribution the energy was 
not appreciably higher, but then no long-range order is 
possible and hence this form can be excluded. In a re- 
cent unbiased projector QMC calculation, the probability 
distribution P(x, y) of the bonds was calculated^ The 
form P(r) ~ r~ 3 was found (with no discernible angu- 
lar dependence). Without a hard-core constraint for the 
VB dimers, the probabilities would clearly be propor- 
tional to the amplitudes; P(x, y) cx h(x, y), and even with 
the hard-core constraint one would expect the two to be 
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ln(L) 

FIG. 3: Log-log plot of the amplitude h(L/2, L/2 — 1) versus 
the system size. Statistical errors are of the order of the size 
of the circles. The line shows the power-law h ~ L~ z . 



strongl y re lated to each other. In fact, as was pointed out 
in Ref. |23j, a wave function with h(r) ~ r~ p does result 
in P(r) ~ r~ p . Our variational calculation confirms that 
indeed the fully optimized h(r) ~ r~ 3 , as demonstrated 
in Fig. [3] using the longest bonds, (x,y) — (L/2, L/2 — 1), 
on the periodic lattices. 

Havilio and Auerbach carried out a VB mean-filed cal- 
culation which gave an exponent p « 2.7*^ The statis- 
tical accuracy in Fig. [3] is perhaps not sufficient to defi- 
nitely conclude that p — 3 exactly, or to exclude p = 2.7, 
from this data alone. However, the QMC study of the 
probability distribution P(R) supports p = 3 to signifi- 
cantly higher precision^ Moreover, Beach has recently 
developed a different mean-field theory which predicts 
p = d + 1 for a d-dimensional system^ There is thus 
reason to believe that r~ 3 indeed is the correct form for 
d = 2. 

For the 4x4 lattice we can compare the variational 
wave function with the exact ground state obtained by 
exact diagonalization. This comparison is most easily 
done by transforming the VB state to the S z basis. Tak- 
ing into account the lattice symmetries, there are 822 
m z = states with momentum k = 0, and the matrix can 
easily be diagonalized. We generate the 8! VB states |T4) 
using a permutation scheme and convert each of them 
into 2 8 S^-basis states with weights ± Y[ h(x, y), and use 
these to calculate the overlap with the exact ground state. 
With the amplitudes normalized by 0) = 1, there is 
only one independent amplitude, h(2, 1), to vary for the 
4x4 lattice. In Fig. 2] we show the overlap as a function 
of h(2, 1). We also indicate the value of h(2, 1) obtained 
in the variational QMC calculation — it matches almost 
perfectly that of the maximum overlap. The best over- 
lap is indeed very high; « 0.9998. It would be interesting 
to see how the overlap depends on the system size. For a 
6x6 lattice, the ground state can also be calculated, using 
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FIG. 4: (Color online) Overlap between the exact 4x4 wave 
function and the VB wave function versus the single indepen- 
dent amplitude h(2, 1). The value of h(2, 1) obtained in the 
variational calculation is indicated. 




- yj 


=0.0 


- yj 


=0.1 


_yj 


=0.2 


-yj 


=0.3 


yj 


=0.4 



0.4 



FIG. 5: (Color online) Overlap beween the exact 4x4 ground 
state and the VB wave function for different values of the frus- 
tration J2/J1. The best amplitudes obtained in variational 
Monte Carlo calculations are indicated with the dashed lines. 



the Lanczos method, but the space of valence bond states 
is too large to calculate the overlap exactly (although it 
could in principle be done by stochastic sampling). 

V. FRUSTRATED SYSTEMS 

We have also studied the Hciscnberg hamiltonian in- 
cluding a frustrating interaction^ 

H=J 1 ^S i -S j + J 2 S i- S i' ( 19 ) 

(i,3) «W» 

where (i, j) and denotes nearest and next-nearest 

neighbors, respectively, and Ji, J2 > 0. Also in this case 
there exists, in principle, a positive-definite expansion of 
the ground state in the valence bond basis. This can 
be easily seen because a negative coefficient fk in Eq. @ 
can be made positive simply by reversing the order of the 
indices in one singlet in that particular state. However, 
no practically useful convention for fixing the order is 
known. We here use the same partition of the lattice 
into A and B sublattice sites as in the non-frustrated 
case and the same sign convention ^ for the singlets. 
We only consider the 4x4 lattice, which, as we will 
show, already gives some interesting information on the 
behavior of the simple amplitude-product wave function 
as the frustration ratio J2/J1 is increased. 

In the exact calculation we can study both positive 
and negative values of h(2, 1), but for now we restrict 
the variational calculation to h(2, 1) > 0, in order to 
avoid the Monte Carlo sign problem caused by negative 
amplitudes^ It should be noted, however, that the sign 
problem here is much less severe than in exact QMC 
schemes^ 3 - and hence there is some hope of actually being 
able to consider mixed signs in variational QMC calcula- 
tions in the VB basis 4^ 



In Fig. [5] we plot the dependence on h(2, 1) of the over- 
lap between the VB wave function and the exact ground 
state for several values of J2I J\- The h(2, 1) correspond- 
ing to maximum overlap decreases as the frustration in- 
creases. For J2/J1 — 0.4 the best overlap occurs for 
h(2, 1) < 0. The optimum overlap decreases significantly 
with h(2, 1) for J2/J1 > 0.3, indicating the increasing ef- 
fects of bond correlations not taken into account in the 
product-form of the expansion coefficients. This deterio- 
ration of the wave function may be related to the phase 
transition taking place in this model at J2/J1 ~ 0.4. 39 
Note, however, that even at J2/J1 = 0.4 the overlap re- 
mains as high as ss 0.996. 

There is a point close to J2/J1 = 0.4 where h(2, 1) van- 
ishes and thus the best wave function for the 4x4 lattice 
contains only bonds of length 1. Beyond this coupling 
the optimum wave function requires a negative h(2, 1). 
It has also been noted previously that wave functions in- 
cluding only the shortest bonds give the best description 
of the ground state in a narrow region of high frustra- 
tion in a model containing also a third-nearest neighbor 
interaction 

We also show in Fig. [5] the values of h(2, 1) obtained 
in the variational calculations. Interestingly, these val- 
ues coincides well with the maximum overlaps only when 
the frustration is weak, showing that the best variational 
state in a given class is not always the best in terms of 
the wave function. 



VI. SUMMARY AND CONCLUSIONS 

In conclusion, we have shown that a variational valence 
bond wave function, parametrized in terms of bond am- 
plitude products, gives a very good description over-all of 
the 2D Heisenberg model. Although this has been known 
qualitatively for a long time^ our study shows that the 
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agreement is quantitatively even better than what was 
anticipated in previous studies. The deviation of the 
ground state energy from the exact value is w 0.06% for 
large lattices; almost 50% better than in previous cal- 
culations where a functional form was assumed for the 
amplitudes.— The sublattice magnetization is correct to 
within w 1% (smaller than the true value). We have also 
shown that the amplitudes for bonds of length r decay as 
r~ 3 for large r, which is the same form as for the prob- 
ability distribution calculated previously— It is also in 
agreement with a recently developed mean-field theory. 28 

By exactly diagonalizing the hamiltonian on a 4 x 
4 lattice, we have also studied the frustrated J1-J2 
model. Not surprisingly, we found that the quality of 
the amplitude-product wave function deteriorates when 
the frustration J%j J\ is increased. However, even at 
J2/J1 — 0.4, i.e., close to the phase transition taking 
place in this model,— the overlap is above 0.996. It would 
clearly be interesting to see how well the amplitude- 
product state works for the frustrated model on larger 
lattices. In this regard we note that Capriotti et aZ.— re- 
cently carried out a variational study of an RVB function 
written in terms of fermion operators^ and found that it 
gave the best description of the ground state of the Ji- 
J2 model at large frustration; J2/J1 ~ 0.5. However, the 
overlap is significantly smaller than what we have found 
here for the 4x4 system at the same level of frustration. 

Although the fermionio^ and bosonic descriptions of 
the VB states are formally equivalent, the fcrmionic wave 
function, as it is normally written, does not span the 
full space possible with the bosonic product state. As a 
consequence, the bosonic description we have used here 
can in practice deliver much better variational wave func- 
tions for non-frustrated systems.— The fermionic descrip- 
tion apparently works better for frustrated than non- 
frustrated systems.— However, if the sign is also opti- 
mized for each amplitude in the bosonic product state 
(which is not easy for large highly frustrated systems, 
however, because of Monte Carlo sign problems 4 ^) it is 
clear that these wave functions should be better than 
the fermionic RVB state considered so far.— Therefore, 



the VB wave function we have studied here should, at 
least in principle, give an even better description of the 
ground state of the frustrated model than the fermionic 
RVB state optimized in Ref. UU. Our results for the 4x4 
lattice, along with the results of Ref. [l3|, suggest that the 
QMC sign problem 4 ^ should be small up to J2/J1 ~ 0.4. 
It may thus even be possible to gain insights into the 
quantum phase transition and the controversial stated 
for J%j J\ > 0.4 with this type of variational wave func- 
tion. 

Another interesting question is how bond correlations, 
which are not included in the wave function considered 
here, develop as the phase transition at J2/J1 ~ 0.4 is 
approached. We are currently exploring the inclusion of 
bond-pair correlations to further improve the variational 
wave function for the Heisenberg model as well as more 
complicated spin models. 

The stochastic energy minimization scheme that we 
have introduced here, which requires only the signs of 
the first energy derivatives, may also find applications in 
variational QMC simulations of electronic systems. Re- 
cently proposed efficient optimization scheme d 33 ' 42 need 
the second energy derivatives, and so our scheme requir- 
ing only first derivatives has the potential of significant 
time savings when the number of variational parameters 
is large. Very recently other powerful schemes also re- 
quiring only the first energy derivatives have been de- 
veloped and have been shown to be applicable to wave 
functions with a large number of parameters 4 s - We have 
not yet compared the efficiencies of these different op- 
timization approaches with the stochastic scheme pre- 
sented here. 
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